A general latent normal model for multilevel data with mixtures of response types is extended in the case of ordered responses to deal with variates having a large number of categories and including count data. An example is analysed by using repeated measures data on child growth and adult measures of body mass index and glucose. Applications are described that are concerned with the flexible prediction of adult measurements from collections of growth measurements and for studying the relationship between the number of measurement occasions and growth trajectories.

In studies of childrens’ growth a common objective is to develop procedures that can be used to predict adult conditions from repeated measurements taken during childhood. Early attempts to do this for adult variables such as stature essentially consider linear prediction models for the response regressed on a set of growth measures taken at a predefined set of ages (usually at whole years of age—see, for example, Tanner et al. (2000)). Such a system, however, lacks flexibility since in practice individuals will be measured at a set of ages that do not typically correspond to the predetermined set. A more flexible approach is to model the growth measures as functions of age, jointly with the adult measurements, by using suitable longitudinal data. The resulting model parameters can then be used to provide a suitable prediction model that will apply, in principle, to any set of growth measurements taken on an individual. In the present paper we consider such a model applied to the prediction of adult glucose levels and body mass index BMI (measured in kilograms per square metre), using measures of weight taken in children under the age of 10 years. These predictions, together with suitable prediction intervals, can then be used as a screening device for the early detection of later problems.

The basic design, therefore, is based on a two-level repeated measures model as described, for example, in Goldstein (2003), chapter 5. Level 1 is defined by the growth measurement occasions and these are clustered within each individual, which is the level 2 unit. The idea is that the growth measurement parameters describing average growth, velocity of growth etc. vary at the individual level and covary among themselves and with the adult outcomes. If we assume multivariate normality, the parameters of the fitted model, notably the ele ments of the individual level covariance matrix, then provide sufficient statistics for the linear prediction of any of the adult measurements from any combination of growth measurements.

One of the problems in repeated measures growth studies is that attrition may be non-random, leading to parameter estimation bias. A further aim of the present paper is to study whether this might be so. A standard procedure for dealing with informative dropout is via mixture models (Hogan and Laird, 1997) where the response vector, i.e. the set of growth measures, and the dropout distribution are jointly modelled. In the present paper we study the issue of informative dropout by modelling the actual number of occasions on which each individual was measured, jointly with the growth parameters and with adult measurements, providing estimates of the relevant correlations.

In the next section we describe a general procedure for these joint models and show how it can be specialized to handle our data.

Goldstein et al. (2009) present a general model (the GCKL model) for multivariate responses that can occur at several levels of a data hierarchy and where the responses can be continuous or discrete. The basic methodology is described in a freely downloadable document from the Centre for Multilevel Modelling Web site (Goldstein et al., 2007): http://www.cmm.bristol.ac.uk.

This work includes, as special cases, existing literature on the joint modelling of discrete and continuous responses, and this literature is referenced there. The GCKL model extends previous work by allowing for a multilevel structure for the data and can also jointly model ordered, unordered and continuous data. The ability to model ordered responses within the GCKL framework provides a very general procedure for any kind of discrete data where we can treat the categories as an ordering, including count data. In the present paper we explore the case where we have a mixture of normal and count data at two levels of a data hierarchy. If we treat count data as an ordered classification Goldstein et al. (2009) provide a procedure to obtain parameter estimates, but this may be inefficient and even numerically unstable, when the number of categories is very large, since it then involves the estimation of a large number of ‘threshold’ parameters. We therefore propose two extensions to the GCKL model: one involves using smoothing functions for the large number of threshold parameters that are involved in the algorithms that were proposed by Goldstein et al. (2009) and the other is based on a latent normal characterization for the Poisson distribution that can model count data that can be assumed to follow such a distribution. Since the growth measure in our data is restricted to weight, it may in practice have limited utility as a sensitive prediction model, but the methodology that we propose will have general applicability.

We now briefly review the two-level GCKL model as it applies to normal and ordered responses. Full details are given in Goldstein et al. (2009). We begin by describing the multivariate normal model and its estimation and in the following section we show how ordered responses can be dealt with by using a transformation to underlying normality.

Let j = 1,…,J index level 2 units and i = 1,…,Ij index level 1 units, nested within the level 2 units. The underlying multivariate normal model is

The superscripts denote the level at which a variable is measured or defined. Here γij(1) is a p1 row vector containing the normal responses that are defined at level 1 for level 1 unit (observation) i nested in level 2 unit j. Without loss of generality, we assume the same set of predictors for each of these level 1 responses. Let X1ij be an f1 row vector that contains observations on these predictor variables for observation i nested in higher level unit j and β(1) is an f1×p1 matrix containing their fixed coefficients. The columns of β(1) correspond to coefficients for each level 1 response variable. Similarly Z1ij is a q1 row vector that contains for observation i nested in higher level unit j values of q1 predictor variables related to the q1 random effects for the level 1 responses. Quite often they may be a subset of the X1ij. The γj(1) is a q1×p1 matrix containing the random effects at level 2 for the level 1 responses, with each column relating to a particular response variable. The random effects may be correlated with each other both within and across level 1 responses. Each response has an associated level 1 residual and these are contained in the p1 row vector eij(1)⁠. Their covariance matrix is Ω1.

Correspondingly, in the second part of the expression yj(2) is a p2 row vector containing the remaining responses that are defined at level 2 and X2j is an f2 row vector that contains predictor variables for higher level unit j, again assumed to be the same for each level 2 response variable. The fixed coefficients of these variables are contained in the f2×p2 matrix β(2). The Z2j is a q2 row vector related to the level 2 random effects for the level 2 responses. The yj(2) is a q2×p2 matrix of correlated level 2 residuals which are also correlated with the level 2 residuals for the level 1 response. The full covariance matrix of all the level 2 residuals is Ω2.

In the examples of this paper we assume that q2 = 1. We shall also assume that the level 1 responses only have simple intercept random effects.

Estimates for the parameters of this model are obtained by using Markov chain Monte Carlo sampling and details are given in Goldstein et al. (2009). Assuming starting values and a burn-in period, given current values, the sampling steps are outlined as follows.

Step 1: sample a new set of fixed coefficients at level 1.

Step 2: sample a new set of fixed coefficients at level 2.

Step 3: calculate a new set of level 2 residuals for level 2 responses by subtracting the fixed predictor X2jβ(2) from the associated normal responses.

Step 4: sample a new set of level 2 residuals for level 1 responses.

Step 5: calculate a new set of level 1 residuals by subtracting the prediction using the fixed effects plus the level 2 random effects, x1ijβ(1)+Z1ijuj(1)⁠, from the associated normal responses.

Step 6: sample a new level 1 covariance matrix.

Step 7: sample a new level 2 covariance matrix.

Missing responses are treated by imputing the corresponding random effect conditional on current values for all the fixed predictors and correlated random effects (Goldstein et al., 2009). In the application that we consider there are no missing responses for the level 2 (individual) unit measurements. We also note that, in the framework of a repeated measures model, there is no requirement to have the same number of level 1 (measurement occasion) units for each level 2 unit (Goldstein (2003), chapter 5).

There are particular issues about the choice of default diffuse prior distributions for the covariance matrices. For example, Goldstein et al. (2009) describe the choice of independent uniform priors for the elements of this matrix, and Browne (2006) suggested that, where maximum likelihood estimates for parameters are readily available from a preliminary fit, it may be reasonable to use an inverse Wishart prior with covariance parameter equal to the corresponding maximum likelihood (or quasi-likelihood) estimate and with degrees of freedom equal to the order of the covariance matrix. We could also use an ‘adaptive prior’ whereby a set of iterations during a secondary burn-in period is averaged and this is used as the prior for a further burn-in period and the final set of iterations. In the examples of the present paper, because we shall be dealing with constrained covariance matrices, we shall need to update the elements of these matrices one by one and this leads us to use a set of independent uniform priors for these elements.

The GCKL model also allows for non-normal responses by formulating a latent normal model where an extra step is inserted to sample an underlying (latent) normal variable from the observed response. Goldstein et al. (2009) show how this can be done for ordered and unordered categorical data. They additionally allow for non-normal continuous data by using a Box–Cox transformation. When such a latent normal variable is sampled, in addition to conditioning on the current values of the parameters and the observations it is also necessary to condition on the correlated responses and this ensures that the latent variables and observed normal variables have a joint multivariate normal distribution at both level 1 and level 2. In the next two sections we outline the procedure for general ordered data and then describe the extension to Poisson data. In all these cases, for ease of exposition, we consider the single-level case. The extension to two levels follows the description that was given above.

Suppose that we have an ordered p-category response indexed by numbers g = 1,…,p. For simplicity we shall use a single-level index i and denote by πgi the probability that observation i occurs in category g. The cumulative distribution is defined as

We consider the probit link model for these as

where ϕ(t) is the density function of the standard normal distribution and (Xβ)i is the value for i of a linear predictor expressed in matrix terms. The αh are known as threshold parameters. In general the linear predictor (Xβ)i will contain higher level random effects and also the terms that arise from conditioning on all correlated random effects from the remaining responses. Since the joint distribution is multivariate normal, as pointed out above, this conditioning is linear and is readily derived from the current value of the covariance matrix. Details are given in Goldstein et al. (2009).

We also note that in equation (2) the linear predictor remains constant over the categories that are defined by the thresholds. Interactions between the threshold parameters and the predictor variables are possible and can be included in straightforward ways.

We can also motivate the latent normal model (2) as follows. If yi is the observed category then yi ⩽ h if and only if a latent variable yi*≤0 where yi*=(Xβ)i−αh+ri,ei~N(0,1)⁠. Thus, we have

If an intercept term is incorporated in the fixed part predictor then there must be a linear constraint on the αh and it is convenient to take α1 = 0 for this.

We can convert the categorical response to the latent normal model by sampling as follows.

For a category 1 response we sample randomly from the standard normal distribution restricted to [−∞,−(Xβ)i].

For a category p response we sample from the standard normal distribution restricted to [αp−1−(Xβ)i,∞].

For every other category h we sample from the standard normal distribution restricted to [αh−1−(Xβ)i,αh−(Xβ)i].

Sampling for the threshold parameters αh is most efficiently carried out by using a Metropolis step (Goldstein et al., 2009).

As mentioned in Section 1, for many kinds of count data, if we treat every observed count value as a separate category, we may have a very large number of threshold parameters to estimate, and furthermore some of the categories may be small and these can cause computational problems. One solution would be to merge adjacent categories where there are small numbers. This, however, suffers from several drawbacks.

First, if there are categories with few observations the associated threshold parameters will not be well estimated and the computational time will be increased. Secondly, although the merging of adjacent categories will not change the parameters being estimated in an additive linear model, this is not so if there are interactions between the threshold parameters and the linear model coefficients. If such interactions exist and are ignored then in general we shall have different values for these coefficients depending on which categories are merged. Likewise, if an interactive model is fitted, and categories involved in these interactions are merged, then different values generally will be estimated. Thirdly, if adjacent categories are merged and this is felt to be justified in terms of the underlying model, it will not allow us to make inferences about the proportions that are associated with the original, unmerged, categories. This restriction does not apply to the smoothing approach that is now described.

In this approach we seek a functional form to describe the threshold values based on a few parameters and the following describes one such approach.

In equation (2) write

where f can be any suitable function. We shall illustrate by using a simple polynomial as a regression function of the category sequence number h typically centred at the mean value. Other choices are possible, including fractional polynomials and regression splines, and we could also use, for example, a set of starting values computed for the thresholds that are derived directly from assigning normal equivalent deviates to the cumulative category probabilities, but we have not investigated these possibilities. A further possibility that is currently being explored (Goldstein, 2008) is a function that consists of a sum of exponential terms that automatically satisfies the monotonicity requirement (see below). We have, omitting subscripts,

We sample each δ in turn by using a Metropolis step with a normal proposal. Using equation (4), with the αh so determined, we then can compute the set {αh} and thus the likelihoods for acceptance or rejection as in the standard algorithm for ordered data (Goldstein et al., 2009). The function f must be strictly increasing over the range of h and this is easily verified at each step: if the order constraints for the set of thresholds are not obeyed we do not update the value. We then sample from the latent normal distribution as in the ordered case as detailed above.

To obtain starting values we first obtain starting values for each separate α as described in Goldstein et al. (2009), and then we fit the regression

to determine the order of polynomial to be fitted and starting values for coefficients. The estimated standard errors for the coefficients can be used for the normal proposal distributions.

Finally, it is possible to have a partially smoothed model where the smoothing takes place only over a contiguous subset of the threshold parameters, e.g. those where numbers are small.

An alternative to the previous section’s treatment of count data is to assume a particular parametric distribution for the counts and we shall consider the common one, namely that the counts follow a Poisson distribution. If a Poisson distribution can be supported by the data then it will be efficient to model the data explicitly assuming such a distribution, since we shall only need to estimate a single parameter. As with general ordered data we can then incorporate this distribution within the GCKL framework that underpins the procedures which are developed in this paper. To do this we need to find a Markov chain Monte Carlo step to sample a latent normal variable and we now show how this can be done.

Consider the Poisson density function

for the first p categories that are observed.

As in equation (2) we write the cumulative distribution

We choose a reference value h = h0 such that

where we use i to index the units.

We note that for a given h0 this is a flexible model relating the Poisson parameter θi for an individual unit to the covariate values. Thus, given (Xβ)i, θi is determined. It is convenient to choose h0 = 0 so that

We sample from the underlying standard normal distribution interval as follows for an observed count h*, where the interval is indicated by brackets {}.

For an observed count in the first category we sample a value from the lower tail of the normal distribution

and for the remaining categories we sample from

The resulting sampled values, just as in the ordered probit model, are thus sampled from the specified standard normal distribution. Since this a novel characterization of the Poisson distribution there is no standard term for our procedure, and we shall refer to it simply as a Poisson latent normal transformation.

To impute any missing values on the original scale we use the inverse transformation (Goldstein et al., 2009).

We note that our models for ordered categories consider the latent normal distribution to be defined by the observed categories, whereas the Poisson model is defined on (0,∞) and so assumes that the distribution exists beyond the categories observed. This suggests that we would expect to obtain different estimates. Since our data counts have a specified minimum count of 1 and a specified maximum possible count, the Poisson assumption may not be reasonable, and this is confirmed in our example. In the Poisson case, in some circumstances, we might consider fitting a truncated distribution that covers only the observed categories, or a particular subset of categories specified in advance. One example is the zero-truncated Poisson distribution where a zero category exists but cannot be observed, but this does not seem to be a reasonable assumption for our data.

Our application uses data that consist of 1000 subjects with serial weight measurements between birth and age 10 years together with adult body mass index BMI and plasma glucose level (millimoles per litre) measured on individuals at around the age of 30 years. During infancy (0–2 years) the maximum number of measurement occasions observed is 9 and during childhood (2–10 years) it is 14. There are 7459 repeated measurements over childhood and 1000 adult BMI and glucose measurements. The data set was made publicly available at the Third International Congress on Developmental Origins of Health and Disease (Gillman et al., 2007).

As explained in Section 1, the joint modelling of the growth data and the adult measurements provides a flexible prediction system. This is the part of the model that is of principal interest. We also include the counts of the number of measurements as a further two, level 2, responses in infancy and in childhood. The principal purpose of this is to see whether the number of measurement occasions is related to characteristics of growth. Thus, we might suppose that children with atypical characteristics may also have atypical growth patterns and that this may be associated with the propensity to continue participation in the study. A strong relationship between the count of the number of measurements and such growth parameters may thus be indicative of non-random dropout. Strictly speaking, a smaller number of measurement occasions is not always indicative of dropout, although it is associated with it in the present data set. A more sensitive measure would be to count the number of measurements within more narrowly defined age groups, but this has not been done for the present data.

The responses in our model for prediction are as follows. At level 1 we have repeated measures of the log(weight) of children up to 10 years (y1ij) with birth weight (x) in kilograms as a covariate. Growth is modelled by using a basic cubic polynomial in age. At adulthood we have measures on the same individuals of log(glucose level) (y2j) and log(BMI) (y3j). We write the model as

Here, in terms of our general notation in equation (1), we have two level 2 response variables yj(2)=(y2j,y3j) and a single level 1 response variable y1ij(1)=y1ij⁠. The term δij, in the child growth model, is an indicator variable taking the value 1 if the weight measurement is taken in infancy (0–2 years) and 0 if taken in childhood (2–10 years) so that the first line of model (10) fits two disjoint polynomial models: one for infancy and one for childhood. The two are linked through the level 2 random effects. This separation was specified by the original investigators who wished to consider the two age ranges separately. An alternative would have been to fit, for example, a regression spline to allow a smooth relationship across the age range, but we have not pursued this here. Age (t) in infancy is centred at 1.0 years and in childhood at 5.5 years. Note that Ωu is of order 6 since the quadratic and cubic terms in model (10) are assumed not to have random coefficients.

The count data are introduced as follows. We have two further responses at level 2, namely a count of the number of measurement occasions in infancy (y4j) and in childhood (y5j). The random (intercept) effects for these (on the latent normal scale) at level 2 are respectively u8j and u9j. The only covariates are the intercepts for these responses.

For the Poisson model we have

where

and

where

where a and b are as defined in expressions (8) and (9).

For the ordered model with threshold parameters we have for infancy

and for childhood

where p4 and p5 are the numbers of categories respectively for the numbers of occasions in infancy and childhood. For the smoothing model, equations (11) and (12) are modified accordingly. Thus the level 2 covariance matrix Ωu has elements σu82=σu92=1 corresponding to the latent normal variables for the ordered responses. The Markov chain Monte Carlo step that samples the remaining elements of Ωu first samples the variance terms followed by the correlations, which are then transformed to covariances.

Table 1 shows the estimates for the fixed parameters for three models. Model A treats the number of measurement occasions in infancy and childhood as ordered categories where each threshold parameter is estimated. Model B smooths the threshold parameters during infancy by using a third-order polynomial and the threshold parameters during childhood by using a fourth-order polynomial. Fitting a higher order polynomial in both cases does not change the other parameter estimates and the higher order polynomial coefficient estimates are smaller than their estimated standard errors. Model C fits a Poisson model. We comment further on these results below.

Fixed effect and threshold estimates for the model defined by expressions (10)–(12)†

Burn-in, 500; sample size, 5000; number of level 1 units, 7459; number of level 2 units, 1000. Standard errors are given in parentheses. For the smoothed model the α-estimates are those derived from the smoothing polynomial at each iteration and the polynomial terms are centred on the mean category (numbered 1,…,8 and 1,…,13 for infancy and childhood respectively). Model A treats the number of measurement occasions in infancy and childhood as ordered categories where each threshold parameter is estimated. Model B smooths the threshold parameters by using a third-order polynomial for infancy and a fourth-order polynomial for childhood, and model C fits a Poisson model.

Fixed effect and threshold estimates for the model defined by expressions (10)–(12)†

Burn-in, 500; sample size, 5000; number of level 1 units, 7459; number of level 2 units, 1000. Standard errors are given in parentheses. For the smoothed model the α-estimates are those derived from the smoothing polynomial at each iteration and the polynomial terms are centred on the mean category (numbered 1,…,8 and 1,…,13 for infancy and childhood respectively). Model A treats the number of measurement occasions in infancy and childhood as ordered categories where each threshold parameter is estimated. Model B smooths the threshold parameters by using a third-order polynomial for infancy and a fourth-order polynomial for childhood, and model C fits a Poisson model.

The chains for the threshold parameters were the least satisfactory in terms of how well they mixed. The burn-in period was chosen after observing the behaviour of these chains, especially for the extreme thresholds, using trial runs with different subsets of parameters. A value of 500 was found to provide a reasonable period for the chains to become stationary. Likewise, we found that 5000 iterations provided reasonable mixing and summary statistics in terms of means and standard deviations that were stable and did not change appreciably with further iterations, using similar trial runs. Further support for the values that were chosen is provided by the fact that models A and B, using different parameterizations, do provide similar estimates for the threshold parameters, albeit with different efficiencies (see below).

The level 2 covariance and correlation matrices are displayed in Table 2 where the variances are on the diagonal, and correlations below the diagonal. It will be seen that, whereas there are moderate correlations between the growth random effects, the correlations of growth random effects with adult measures are not large, so the usefulness of weight measures for prediction of the adult measures is limited. To illustrate the prediction procedure we consider predicting log(BMI) by using two childhood weight measurements, together with birth weight.

Level 2 variances and correlations for model (10)–(12)†

The response variables are in the following order: overall intercept; infancy indicator; infancy slope; childhood slope; log(glucose); log(BMI); infancy count; childhood count.

Level 2 variances and correlations for model (10)–(12)†

The response variables are in the following order: overall intercept; infancy indicator; infancy slope; childhood slope; log(glucose); log(BMI); infancy count; childhood count.

We write the ‘raw’ residuals for the two weight measurements at times t1j and t2j and the adult log(BMI) from model (10) as

This defines a three-variate normal distribution that is given by

From this normal distribution we can derive the expected value of any of the responses given the other two, which is just a linear (regression) predictor, which for z3j has the form

And the regression coefficients are derived from the above expression, namely

where sample estimates, such as those given in Table 2, are substituted. The estimated covariance matrix of the coefficients is obtained from the inverse matrix given in the expression above and this then enables us to place confidence intervals around the predicted values Ζ^3j⁠. The predicted value of log(BMI) is given by combining Ζ^3j with Ζ^1 so that the predicted value of BMI becomes

The count responses do not follow a Poisson distribution very closely, as is shown in Fig. 1; there is a non-significant χ2 goodness-of-fit statistic for infancy with a significant value for childhood. This is also reflected in some of the estimates in Table 1 (see below) and the differences in the correlations in the last two rows of the correlation matrices in Table 2.

Observed probability distributions for infancy and childhood number of measurement occasions (categories coded starting from 1; Poisson models were fitted): (a) infancy (χ2 = 6.85 (8 degrees of freedom), P>0.05); (b) childhood (χ2 = 39.1 (13 degrees of freedom), P

Google Scholar

WorldCat

Goldstein, H. (

Google Scholar

Google Preview

WorldCat

Goldstein, H. (

Goldstein, H., Carpenter, J., Kenward, M. and Levin, K. (

Google Scholar

WorldCat

Goldstein, H., Steele, F., Rasbash, J. and Charlton, C. (

Google Scholar

Google Preview

WorldCat

Hogan, J. W. and Laird N. M. (

Google Scholar

WorldCat

Google Scholar

Google Preview

WorldCat

Tanner, J. M., Healy, M. J. R., Goldstein, H. and Cameron, N. (

Google Scholar

Google Preview

WorldCat